Polaron effective mass from Monte Carlo simulations 
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A new Monte Carlo algorithm for calculating polaron effective mass is proposed. It is based on 
the path-integral representation of a partial partition function with fixed total quasi-momentum. 
Phonon degrees of freedom are integrated out analytically resulting in a single-electron system with 
retarded self-interaction and open boundary conditions in imaginary time. The effective mass is 
inversely proportional to the covariance of total energy calculated on an electron trajectory and the 
square distance between ends of the trajectory. The method has no limitations on values of model 
parameters and on the size and dimensionality of the system although large statistics is required 
for stable numerical results. The method is tested on the one- dimensional Holstein model for which 
simulation results are presented. 

PACS numbers: 63.20.Kr, 71.38.+i 



The polaron problem is a classic problem of condensed- 
matter physics pl|| . In recent years it has attracted ad- 
ditional attention in relation to high-temperature super- 
conductivity ||. Currently active theoretical research 
is being conducted by various methods on models with 
strong electron-phonon phonon interaction as well as on 
models with both electron-phonon and electron-electron 
interaction such as Hubbard- Holstein and t — J- Holstein 
models. The Monte Carlo simulation method is widely 
used in these studies [BJ because it provides unbiased in- 
formation about thermodynamic and dynamic properties 
of the system even in regions of model parameters hardly 
accessible by standard theoretical tools. 

The polaron problem was among the first applications 
of the Quantum Monte Carlo method after its introduc- 
tion in condensed-matter physics in the Seventies. In 
simulations of many-polaron systems either both the po- 
sitions of the electrons and the ions are simulated PJ3] 
or the electrons are integrated out analytically and only 
the phonon subsystem is simulated B. For the single- 
polaron problem, De Raedt and Lagendijk || developed 
an effective algorithm in which phonon degrees of free- 
dom are integrated out analytically by the Feynman tech- 
nique [J9|. This results in a single-electron problem with 
retarded self-interaction along the imaginary-time direc- 
tion. Since a system with only one degree of freedom can 
be easily simulated by the Metropolis Monte Carlo algo- 
rithm fllOf this method provides accurate results for ther- 
modynamic polaron properties such as internal energy 
and specific heat and on other static quantities of in- 
terest, for instance the electron-phonon correlation func- 
tion. Unfortunately, the estimation of dynamic proper- 
ties with this method, such as the polaron spectrum, is 
much more difficult, the reason being that this would re- 
quire an analytic continuation of simulation results from 
the imaginary-time to the real-time domain. This pro- 
cedure is mathematically ill-posed and is thus very sen- 
sitive to the statistical noise in the input data. In pre- 



vious work |11[| we have applied the singular-value de- 
composition technique to carry out such a continuation. 
However, due to the narrowing effect jllEJ the polaron 
band is expected to be very narrow, especially in the 
strong-coupling regime, and the bandwidth smaller than 
the accuracy with which the polaron spectrum could be 
obtained by analytic continuation. 

In this paper we show that despite these difficulties 
the polaron effective mass, which is defined as an in- 
verse second derivative of internal energy with respect to 
quasi-momentum in the limit of zero temperature, can be 
obtained directly from Monte Carlo simulations without 
going to real time. This possibility is based on a spe- 
cial kind of fluctuation-dissipation relation and resembles 
Feynman's calculation of the Frohlich polaron effective 
mass M and the method of calculation of the superfluid 
fraction in liquid helium fll2| ]. 

Our starting point is the fact that total quasi-momen- 
tum P is a conserved quantity in an electron-phonon sys- 
tem. Since states with different P do not mix during the 
evolution it is meaningful to calculate a partial partition 
function Zp which is a sum over states with fixed P 
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Here H is the hamiltonian, (3 — 1/fesT, n numbers states 
of a complete basis in the Hilbert space and P n is the 
total quasi-momentum of state \n). The total partition 
function is obviously Z — J2p Z-p . The total Hilbert 
space is a direct product of the one-electron {k} and 
the phonon Hilbert spaces. The latter in turn can be di- 
vided into subspaces {w^} with fixed total phonon quasi- 
momentum Q. We rewrite (H) as follows 
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We have used a momentum-space basis because the con- 
servation of quasi-momentum is conveniently written in 



momentum space. However, matrix elements of the 
evolution operator exp(— (3H) are easily calculated in 
real space jq]. Therefore, our first goal is to transform 
Eq.(||) to a real-space basis. We do this in two steps. 
Firstly, introducing the Wannier basis for the electron 
r) = TV -1 / 2 ^ k e lkr |k) and exercising the delta- function 
one gets 
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where N is the total number of sites in the lattice. Sec- 
ondly, we introduce a new basis for the phonon states 
in terms of ion configurations |{£i}}, i = 1,...,N. Such 
a configuration represents a state where i-th ion is dis- 
placed from its equilibrium position by distance £». (For 
simplicity we assume only one phonon degree of freedom 
per lattice site). These states are normalised as follows 
({&}|{J» = Ilti^te - #) and resolve the identity: 
I = /_ D£, \{£i}){{Ci}\- Insertion of the identity twice 
into Zp yields 

1 f°° 

Zp = ^E elP(r ~'l £W<r',{e0|e-^|{eO,r) • W 

r,r' J -°° 

W = J2e- iQ(r '- r) Y,({&\ n %)( n Ph\{®)- (3) 

Q rfl h 

The quantity W can now be simplified by the follow- 
ing argument. Let us expand configuration |{£i}} in 
a Fourier series over states with definite wave vectors 

q: |{6» = ^- 1/2 E q a q e iqRi |q>- By |q) we mean a 
state in which ion displacements are arranged in a plane 
wave with wave vector q and the unit amplitude. Now 
comes an important observation which is crucial for the 
whole method. The only component of the latter expan- 
sion which survives projection onto (nr h \ is the one with 
q = Q. Although |q) is a wave of classical displace- 
ments with no definite phonon occupation numbers its 
expansion in quantum states \n p h) would contain only 
those with total quasi-momentum Q = q (this is sim- 
ply because both the classical and the quantum states 
must belong to the same irreducible representation of 
the translation group). It is convenient to take care of 
this property by introducing an extra summation over 
the lattice <5 Q , q = iV" 1 £m e4m(Q ~ q) - Now configura- 
tions ({£i} and {£,'}} are automatically projected on 
the subspace with total quasi-momentum Q and the sum 
E n « I'ViK'V/J in Ec f(|) ma y be replaced by the unity 

ph " " 

operator. With these transformations W takes the form 
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Deriving this result we have first summed over Q, then 
transformed from |q) states back to |{£j}) states and fi- 
nally used normalization of the latter. Substitution of 
Eq.(Q) into Eq.(g) and integration over {^} leads now to 
the final result 

1 r°° 

Zp = ^E eiP(r ~ r) / £>£<r',{6- r ' +r }|e-^|{6},r>. 

r,r' J ~°° 

(5) 

Thus, unlike the total partition function, the par- 
tial partition function with fixed P is a sum over all 
imaginary-time trajectories with open boundary condi- 
tions. Moreover, the boundary conditions for the elec- 
tron and the phonons are correlated. If the final position 
of the electron (i.e. at imaginary time r = 0) is shifted 
from the initial one (at r = (3) by Ar = r' — r then the 
final phonon configuration should be obtained from the 
initial one by shifting the latter by the same Ar along 
the lattice. Each Zp receives contributions from trajec- 
tories with all possible Ar. Obviously, Eq.(g) satisfies the 
condition ^ p Zp = Z . 

The importance of the conservation of total quasi- 
momentum in the polaron problem was first recognized 
by Lee, Low and Pines [[13[ who used a canonical transfor- 
mation to eliminate electron coordinates from the prob- 
lem. Here we follow a different strategy, namely to elim- 
inate phonon degrees of freedom in the spirit of the De 
Raedt and Lagendijk's method. The resulting single- 
electron problem can be simulated by Monte Carlo. Un- 
fortunately, the complex factor e lP ^ r _r ' appearing in 
Eq.(H) makes it impossible to simulate the system at ar- 
bitrary P due to the complex trajectory weight. How- 
ever, there is no sign-problem at P = 0. This fact may 
be used, for instance, for more accurate estimation of po- 
laron ground state energy from the Monte Carlo method. 

In this paper we would like to point out the possibility 
of calculating polaron effective mass, which arises from 
equation (g). One can define the partial internal energy 
in the usual manner as Up = —Zp dZp/d/3. In the 
limit of zero temperature Up is reduced to the lowest 
eigenvalue with quasi-momentum P. Then the inverse 
effective mass is (up to factor % 2 ) a second derivative of 
Up with respect to a chosen component of P: 1/rn* = 
d 2 Up/dPl\p =0 . From Eq.(g) we have 

-L ex - (((r' a - r a fE) - ((r' a r a f) (E)) | P=0 (6) 



where E is the value of the energy estimator calculated on 
the trajectory which begins at r and ends at r' and (...) 
means the Monte Carlo average over an ensemble of such 
trajectories, e.g. (E)\p = Up. Thus, the polaron effective 
mass is inversely proportional to the covariance of energy 
of the trajectory and the square distance between ends of 
the trajectory. Its value can be obtained directly from 
Monte Carlo simulations of the equilibrium system with 
open boundary conditions at P = 0. 

We point out that, on the semi-intuitive level, the re- 
lation between the width of the electron trajectory and 
the polaron effective mass was used previously to demon- 
strate the increase of the latter ||. In the present paper, 
however, such a relation, Eq.(|6|), is derived rigorously. 

In order to test this method we now consider the one- 
dimensional Holstcin polaron problem [|l4| . In this model 
the phonon subsystem is a chain of non-interacting oscil- 
lators with frequency lu and the electron interacts linearly 
only with the oscillator it currently occupies. The model 
hamiltonian reads 

H = H 1 +H 2 + H 3 (7) 

Hi = -t £ (4 + ^ + 4«+i) ; ^ = £ H 

i i 

#3 = E^r^ ~-9E c l c ^- 



Here H\ represents electron kinetic energy, t being the 
nearest-neighbour hopping amplitude, m is the oscillator 
reduced mass, £j is the displacement of i-th oscillator, 
pi = —ihd/d^i and g in the electron-phonon coupling 
constant. The rest of the notation is standard. In what 
follows periodic boundary conditions in real space are 
assumed. 

To be able to use the Monte Carlo method one needs 
a path-integral representation of the matrix element G 
which appears in Eq.(pJ). Following the standard proce- 
dure p|q| we divide the imaginary-time dimension (whose 
extent is /3) into M 3> 1 intervals, insert the resolution of 
the identity M — 1 times, then use the Trotter decomposi- 
tion and evaluate all the matrix elements of the operator 
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I(x i+ i -Xj) = — Ecos [k n (x j+1 - Xj )] e 2TCOsk ". (8) 
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Here {xj} is an electron trajectory in imaginary time 
with boundary conditions Xq — x, Xm = x 1 ] k n are single- 
particle momenta allowed in a chain of N sites: k n = 
2im/N , n = 0, . . . N—l; r = (3t/M and S p h is the phonon 
action 
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In Ref. p| polaron thermodynamics were studied. There 
the following result was obtained for the integral over £s 
with periodic boundary conditions ^m = £jo 
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where Z p h is the partition function of free phonons and 

F(j — j') is the memory function 
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with g 2 = h g 2 jmt? and ui = Tuo/t. In our case one has 
to perform the integration over £s with twisted bound- 
ary conditions £i+ x '—x,M = £i,o- Having done this, how- 
ever, we found that, as long as the condition finco S> 1 
is satisfied, the final result for i' - i / acquires only 
exponentially small corrections relative to Eqs.([10|)-([ll|). 
The details of the calculation are cumbersome and will 
be published elsewhere p5[ . 

Thus, in the low-temperature limit the path-integral 
representation of the partial partition function for the 
model (R) has the form 
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Introducing the definition (A) = Zq T^s x .\Ap{{%j}) 
the expression for the effective mass can be brought to 
the form 

r ^ = -\{{{x'-xfE)-{{x'-xf){E)) (13) 
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where mo = h 2 /2ta 2 is the bare electron mass (a is the 
lattice spacing). 

The inverse effective mass was calculated by the 
Metropolis Monte Carlo algorithm for a chain of N — 
1024 oscillators. The frequency was Co = 1 and the in- 
verse temperature was taken f3 = 15 t~ l to make sure the 
low-temperature limit is satisfied and Eq.(n2|) is valid. 
The number of imaginary-time slices was varied from 
M = 120 to M = 180 to detect a possible M-dependence 



of the simulation results. We observed such a dependence 
in the small-polaron regime g > 2.5. For these values of 
g a l/Af 2 -scaling was used to extrapolate to M = oo. 
In the beginning of each series of measurements we used 
50, 000 • M single-particle steps to warm up the system. 
After that consecutive measurements were taken every 
M steps. For each set of parameters from six to nine 
series of 2 • 10 6 and 3 • 10 6 measurements were made. 




FIG. 1. Inverse polaron effective mass vs coupling constant 
for the one-dimensional Holstein model. Phonon frequency 
u> = 1 and the inverse temperature is /3 — 15 £ -1 . The dashed 
line is mo/m* = exp(— g 2 /2Cj 3 ). 

Results of the simulations averaged over the differ- 
ent series are shown in Fig.l. The effective mass in- 
creases exponentially with the coupling constant, reach- 
ing m* ~ 100 mo in the small polaron regime g > 2.5. 
In spite of the large number of measurements statistical 
errors are not small but decrease with g. This reflects 
the nature of the object simulated. In the absence of 
electron-phonon coupling the electron trajectory is very 
flexible and its ends fluctuate strongly. As the interaction 
is turned on the trajectory becomes more and more rigid. 
This is because trajectories with straight segments ac- 
quire exponentially large weights (see Eq.p3)) and dom- 
inate the sampling. Also shown in Fig.l is the result of 
Ref. [M ( see a l so 050)0) f° r t ne inverse effective mass 
(or the polaron bandwidth) too/™* = exp(— g 2 /2a> 3 ). 
In the intermediate-coupling regime (1.0 < g < 2.5) at 
u> = 1.0 the Monte-Carlo polaron mass is lighter than 
exp(g 2 /2), which is in agreement with exact diagonaliza- 
tion of finite clusters U% . 

In conclusion, we have developed a Monte Carlo algo- 
rithm for calculating polaron effective mass. It is based 
on the representation of the partition function with fixed 
quasi-momentum as a path-integral with open boundary 
conditions in imaginary time. The boundary conditions 
for the electron and phonon subsystem are correlated. 
After analytical elimination of phonon degrees of free- 



dom the resulting single-electron system with retarded 
self-interaction can be simulated by Monte Carlo at zero 
total quasi-momentum. The polaron effective mass turns 
out to be inversely proportional to the covariance of the 
energy of the electron trajectory and the square distance 
between ends of the trajectory. The method does not 
impose any particular limitations on the size and dimen- 
sionality of the system and on values of the model param- 
eters provided a low enough temperature can be reached 
to satisfy the condition f3hu> S> 1. We have tested the 
method on the one-dimensional Holstein model and ob- 
tained physically reasonable results. Statistical errors do 
not appear to be small and large statistics might be re- 
quired to get stable numerical results. 

We are thankful to Prof. A. S. Alexandrov for attract- 
ing our attention to this problem and to V. Kabanov and 
E. Klcpfish for numerous and helpful discussions. 
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